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ABSTRACT 


The dynamics and structure of accretion disks, which accumulate a vertical magnetic field in their centers, are in- 
vestigated using two- and three-dimensional MHD simulations. The central field can be built up to the equipartition 
level, where it disrupts a nearly axisymmetric outer accretion disk inside a magnetospheric radius, forming a mag- 
netically arrested disk (MAD). In the MAD, the mass accretes in the form of irregular dense spiral streams, and the 
vertical field, split into separate bundles, penetrates through the disk plane in low-density magnetic islands. The ac- 
creting mass, when spiraling inward, drags the field and twists it around the axis of rotation, resulting in collimated 
Poynting jets in the polar directions. These jets are powered by the accretion flow with an efficiency of up to ~1.5% 
(in units of Mic). The spiral flow pattern in the MAD is dominated by modes with low azimuthal wavenumbers 
m ~ 1 —5 and can be a source of quasi-periodic oscillations in the outgoing radiation. The formation of the MAD and 
the Poynting jets can naturally explain the observed changes of spectral states in galactic black hole binaries. Our 
study is focused on black hole accretion flows; however, the results can also be applied to accretion disks around 
nonrelativistic objects, such as young stellar objects and stars in binary systems. 


Subject headings: accretion, accretion disks — black hole physics — galaxies: jets — gamma rays: bursts — 
instabilities — ISM: jets and outflows — magnetic fields — MHD — stars: oscillations — 


turbulence 


1. INTRODUCTION 


Accretion disks can carry small- and large-scale magnetic fields. 
The small-scale field (£ < R, where £ is the field scale length and 
R measures the radial distance from the disk center) can be lo- 
cally generated by an MHD dynamo (Brandenburg et al. 1995; 
Stone et al. 1996) supported by the turbulence that results from 
the magnetorotational instability (MRI; Balbus & Hawley 1991). 
This field can provide an outward transport of angular momen- 
tum in the bulk of the disk with the help of local Maxwell stresses 
(Shakura & Sunyaev 1973; Hawley et al. 1996). The large-scale 
field (£ > R) is not likely to be produced in accretion disks (how- 
ever, see Tout & Pringle 1996); rather, it can either be captured 
from the environment and dragged inward by an accretion flow 
(Bisnovatyi-Kogan & Ruzmaikin 1974, 1976) or inherited from 
the past evolution of an accreting object (see § 4). The large-scale 
field can remove the angular momentum from accretion disks 
by global Maxwell stresses through the magnetized disk corona 
(Kónigl 1989). 

A large-scale bipolar field, unlike a small-scale field, cannot 
dissipate locally due to magnetic diffusivity and cannot be absorbed 
by the central black hole. In the case of inefficient outward diffusion 
of the bipolar field through the disk (see Narayan et al. 2003; Spruit 
& Uzdensky 2005; Bisnovatyi-Kogan & Lovelace 2007; and for 
another possibility, see van Ballegooijen 1989; Lubow et al. 1994; 
Lovelace et al. 1994; Heyvaerts et al. 1996; Agapitou & Papaloizou 
1996; Livio et al. 1999), this field is accumulated in the innermost 
region of an accretion disk and forms a “magnetically arrested 
disk," or MAD (Narayan et al. 2003). The MAD consists of two 
parts: the outer, almost axisymmetric, Keplerian accretion disk 
and the inner, magnetically dominated region, in which the accu- 
mulated vertical field disrupts the outer disk at the magnetospheric 
radius R,, ~ 8GM p/B? , where M is the central mass, p is the ac- 
cretion mass density, and B is the magnetic induction. 

It is believed that the large-scale bipolar field in accretion 
disks is responsible for the formation of jets observed in a large 
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variety of astrophysical objects (e.g., Livio et al. 2003). The mag- 
netically driven jets can be of two types, basically depending on 
the mass load of the disk matter (e.g., Lovelace et al. 2005): 
Poynting jets and hydromagnetic jets, which have, respectively, 
small and large mass loads. The hydromagnetic jets can be formed 
by two mechanisms: the magneto-centrifugal mechanism ( Blandford 
& Payne 1982; Kónigl & Pudritz 2000) and the toroidal-field pres- 
sure generated by the disk rotation ( Lynden-Bell 2003; Kato et al. 
2004). The magneto-centrifugal mechanism produces relatively 
wide outflows and, to be consistent with observations of the col- 
limated jets, requires an additional focusing mechanism. Jets 
driven by the toroidal-field pressure can have a high degree of 
collimation, but these Jets are known to be kink unstable ( Eichler 
1993; Appl 1996; Spruit et al. 1997). Poynting jets are naturally 
self-collimated and are expected to be marginally kink stable (Li 
2000; Tomimatsu et al. 2001). These jets can originate in the inner- 
most region of accretion disks and are powered either by the disks 
themselves (Lovelace et al. 1987, 2002) or by rotating black holes 
(Blandford & Znajek 1977; Punsly 2001; see also Takahashi et al. 
1990; Komissarov 2005; Hawley & Krolik 2006; McKinney 2006). 

Our study is based on two- and three-dimensional MHD sim- 
ulations and has two main goals. First, we investigate the dynamics 
and structure of MADs. We show that MADs are formed in ac- 
cretion flows, which carry inward large-scale poloidal magnetic 
fields. Inside the magnetospheric radius Rm, the matter accretes as 
discrete streams and blobs, fighting its way through the strong 
vertical magnetic field fragmented in separate bundles. Because of 
rotation, the streams take spiral shapes. Second, we demonstrate a 
link between the existence of MADs and the production of pow- 
erful Poynting jets. These jets should always be generated in MADs 
because of the interaction of the spiraling accretion flow with the 
vertical magnetic bundles, which, as a result, are twisted around the 
axis of rotation. 

This paper extends the work of Igumenshchev et al. (2003) 
by studying in more detail radiatively inefficient accretion disks 
with poloidal magnetic fields. We employ a new version of our 
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3D MHD code, which can be utilized in multiprocessor simula- 
tions. The paper is organized as follows. We describe the equa- 
tions, the numerical method used, and the initial and boundary 
conditions in 8 2. We present our numerical results in $ 3, and 
discuss and summarize them in 8 4. 


2. NUMERICAL METHOD 


We simulate nonradiative accretion flows around a Schwarzschild 
black hole of mass M using the following equations of ideal MHD: 


dp 


qo (1) 
O82 pk vo +- (Y xB) xB (2) 
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where v is the velocity, P is the gas pressure, ® is the gravitational 
potential, e is the specific internal energy of gas, and q is the total 
energy flux per unit square (see, e.g., Landau & Lifshitz 1987). We 
adopt the ideal gas equation of state, 


P = (y — Upe, (5) 
with an adiabatic index y = 5/3. We neglect self-gravity of the 


gas and employ a pseudo-Newtonian approximation (Paczyński 
& Wiita 1980) for the black hole potential 


(6) 


where R, = 2GM. /c? is the gravitational radius of the black hole. 
No explicit resistivity and viscosity are applied in equations (2)— 
(4). However, because of the use of the total energy equation 
(eq. [3]), the energy released due to the numerical resistivity and 
viscosity is consistently accounted for as heat in our simulations. 

We study the radiatively inefficient regime of accretion, so that 
the energy equation (eq. [3]) does not include any terms represent- 
ing radiation cooling. This regime is appropriate in the two lim- 
iting cases of the very optically thick and very optically thin flows, 
which correspond to the sufficiently sub- and super-Eddington 
mass accretion rates, respectively (see, e.g., Igumenshchev & 
Abramowicz 2000). 

The MHD equations (1)—(4) are solved by employing the time- 
explicit Eulerian finite-difference method, which is an extension 
to MHD of the hydrodynamic piecewise-parabolic method by 
Colella & Woodward (1984). We solve the induction equation 
(eq. [4]) using constrained transport (Evans & Hawley 1988; 
Gardiner & Stone 2005), which preserves the V - B = 0 condi- 
tion. In our method, we employ the approximate MHD Riemann 
solver of Li (2005). Test simulations have shown that this solver 
is robust and provides good material interface tracking. 

We use spherical coordinates (R, 0, à). Our 3D numerical grid 
has 182 x 84 x 240 zones in the radial, polar, and azimuthal direc- 
tions, respectively. The radial zones are spaced logarithmically 
from Ri, = 2R, to Rout = 220R,. Both polar cones with an open- 
ing angle of 7/8 are excluded. Therefore, the polar domain ex- 
tends from 0 = 7/16 to 157/16. The grid resolution in the polar 
direction is gradually changed from a fine resolution around the 
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equatorial plane to a coarse resolution near the poles (with a 
maximum-to-minimum grid size ratio +3). The azimuthal zones 
are uniform and cover the full 27 range in $. The absorption con- 
dition for the mass and the condition of the zero-transverse mag- 
netic field are applied in the inner and outer radial boundaries, 
providing that the mass and field can freely leave the computational 
domain through these boundaries. In the boundaries around the 
excluded polar regions, we apply reflection boundary conditions, 
which means that no streamlines or magnetic lines can go through 
these boundaries. 

Atthe beginning of our simulations, the computational domain 
is filled with a very low-density, nonmagnetized material. The sim- 
ulations are started in 2D, assuming axial symmetry, with an injec- 
tion of mass in a slender torus, which is located in the equatorial 
plane at Rinj = 210R,;, i.e., quite close to Rout. This mass has a 
Keplerian angular momentum and a specific internal energy of 
€inj = 0.045 GM/Riy. After an initial period of simulation with- 
out magnetic fields, the mass forms a steady thick torus, which 
has its inner edge at R ~ 150R, and whose outer half is truncated 
at Rout- The torus contains a constant amount of mass and is in dy- 
namic equilibrium: all the injected mass flows outward through 
Row after circulation inside the torus. No accretion flow is formed 
at this point. This hydrodynamic, steady, thick torus is used as an 
initial configuration in our MHD simulations. 

The MHD simulations are started at t = 0 from the steady, 
thick torus by initiating the injection of a poloidal magnetic field 
into the slender torus at Rinj. The numerical procedure for the 
field injection is similar to that described by Igumenshchev et al. 
(2003) except for one modification: now the strength of the in- 
jected field can be limited by setting the minimum inj, which is 
the ratio of the gas pressure to the magnetic pressure at Rinj. This 
modification allows us to better control the rate of field injection. 
The entire volume of the thick torus is filled by the field during 
about one orbital period, forb, estimated at Rinj. From this moment, 
t ~ torb» the formation of the accretion flow begins as a result of 
the redistribution of the angular momentum in the torus due to 
Maxwell stresses. In the following discussion, we employ dimen- 
sionless time, which is obtained using the normalization timescale 
torb, 1.€., £ — t/torb- 


3. RESULTS 


We present the results of combined axisymmetric 2D and non- 
axisymmetric 3D simulations. The initial evolution in our mod- 
els has been simulated in 2D. This allows us to consider longer 
evolution times than those that can be obtained in 3D simulations, 
because of the larger requirements for computational resources in 
the latter case. We initiate 3D simulations starting from developed 
axisymmetric models. The results have shown that nonaxisym- 
metric motions are not very important in the outer parts of the con- 
structed accretion flows, and therefore the use of 2D simulations 
on the initial evolution stages is a reasonable simplification. 

We consider three models, which differ in their rates of field 
injection as determined by inj = 10, 100, and 1000, and we will 
refer to these as models A, B, and C, respectively. All other prop- 
erties of the models, including the injection radius Rinj, internal 
energy ei, and numerical resolution, are the same (see § 2). 


3.1. Accretion Flows 


The initial axisymmetric development of the models is qual- 
itatively similar: the inner edge of the thick torus is extended 
toward the black hole, forming relatively thin, almost Keplerian 
accretion disks. The time of the disk development is varied, de- 
pending on the strength of the injected field. The accretion of the 
mass into the black hole begins at t ~ 0.7 in model A, ~1.3 in 
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model B, and ~4.2 in model C (time is given in units of the or- 
bital period at Rinj; see 8 2). At this stage, the evolution of the disks 
is governed mainly by global Maxwell stresses produced by the 
poloidal field component. This component is advected inward 
with the accretion flow and, because of the disks' Keplerian rota- 
tion, generates relatively strong toroidal magnetic fields localized 
above and below the midplane. These toroidal fields form a highly 
magnetized disk corona with a typical 9 ~ 0.01. 

Model B, and especially model C, demonstrate the develop- 
ment of 2D MRI. This development is similar to that observed 
by Stone & Pringle (2001), in particular in their “run C." We 
have found the origin of the channel solution (see Hawley & 
Balbus 1992) in the central regions of models B and C. This solu- 
tion consists of oppositely directed radial streams and is the char- 
acteristic feature of the axisymmetric nonlinear MRI (Stone & 
Norman 1994). Analysis of the models shows that the channel 
solution is developed when the wavelength 2 = 2xV4//3Q of 
the fastest growing mode of the MRI is well resolved on the nu- 
merical grid, i.e., 4 Z 5A x, where Ax is the grid size, V4 is the 
Alfvén velocity, and €? is the angular velocity. Model A shows no 
indication of the MRI, which can be attributed to the strong mag- 
netic fields, which suppress the instability. In this model, the es- 
timate of À typically exceeds the disk thickness. Model A has some 
resemblance to the nonturbulent “run F" of Stone & Pringle 
(2001). 

In spite of the mentioned similarities to the results of Stone & 
Pringle (2001), our models show different behavior over long 
evolution times. Our simulation design with a permanent injec- 
tion of mass and magnetic field results in accretion disks, which 
accumulate a poloidal field in the center and form MADs. The 
models of Stone & Pringle (2001; also De Villiers et al. 2003; 
Hirose et al. 2004; McKinney & Gammie 2004; Hawley & Krolik 
2006) did not form MADs and did not show a history of perma- 
nent accretion over a long time, probably because of the initia- 
tion of simulations from static magnetized tori, which contain a 
limited amount of mass and magnetic flux of one sign. 

Our simulation design, in principle, could result in a stationary 
state, in which the magnetospheric radius R,, moves all the way 
outward to the injection radius Rinj. This state, however, is not 
physically relevant because it is determined by an assumed value 
of Rinj and conditions of the mass and field injection. Note that such 
a stationary state was obtained in simulations by Igumenshchev 
et al. (2003) in their model B. In the present simulations, because 
we use different (relatively smaller) field injection rates and fol- 
low relatively shorter evolution periods, we have never reached 
this state. 

In the following discussion, we will concentrate on the results 
describing the formation, evolution, and structure of the inner, 
magnetically dominated regions of the MADs. Other aspects of 
our results, including the structure and dynamics of the outer, 
Keplerian regions of the disks, will be reported elsewhere. 

Figure 1 shows example snapshots of the axisymmetric den- 
sity distribution in model B from 2D simulations at two successive 
moments, 7 = 5.1153 and 5.1458. The accretion flow is nonuni- 
form because of the development of turbulence. The turbulence 
results from the combined effect of the MRI and current sheet 
instability. The latter instability locally releases heat due to recon- 
nections of oppositely directed toroidal magnetic fields. The recon- 
nection heat makes a significant contribution to the local energy 
balance in the central regions of the flow because of the relatively 
high energy density of the field, which is comparable to the grav- 
itational energy density of the accretion mass. The thick disk 
structure observed in Figure 1 is explained by convection motions 
supported by the reconnection heat. Note that the case in which 
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the turbulence is supported by only convection motions from the 
reconnections, without the effects of rotation and MRI, has been 
demonstrated in simulations of spherical magnetized accretion 
flows (Igumenshchev 2006). In the case of disk accretion, almost 
axisymmetric convection motions, similar to those found here, 
have been observed in 3D models with toroidal magnetic fields 
(Igumenshchev et al. 2003). The convection motions in models B 
and C make these models relevant to convection-dominated ac- 
cretion flows (Narayan et al. 2000; Quataert & Gruzinov 2000). 
Our simulations show that the poloidal field is transported inward 
in axisymmetric turbulent flows and accumulated in the vicinity 
of the black holes. When the central poloidal field reaches a 
certain strength (about equipartition with the gravitational en- 
ergy of the accreting mass), the accretion flow becomes unstable 
(Narayan et al. 2003). In axisymmetric simulations, the insta- 
bility takes the form of cycle accretion, in which more extended 
periods of halted accretion (see Fig. 1a) are followed by relatively 
short periods of accretion (see Fig. 15). In the case of the halted 
accretion period, the inner accretion disk is truncated at the mag- 
netospheric radius Rm, which is ~15R, in Figure la. The pressure 
of the strong central vertical field (see Fig. 2a) prevents the mass 
accumulated at R,, from falling into the black hole. The accretion 
begins as soon as the gravity of the accumulated mass overcomes 
the magnetic pressure. During the accretion period, the whole 
magnetic flux, which is localized inside R, in the period of halted 
accretion, is now localized on the black hole horizon (see Fig. 2b). 
Note that similar structural features of the inner MHD flows in 
accretion disks related to the model of gamma-ray bursts were 
discussed by Proga & Zhang (2006). 

Figure 3 illustrates the time dependence of the accretion flow 
in model B, showing the evolution of the mass accretion rate Min 
and magnetic fluxes in the midplane inside the five specific radii: 
210R; (=Rinj), 100R,, 50R,, 25R,, and 2R, (=Rin). This figure 
shows the evolution, which has been simulated in 2D from t = 0 
to 2.14 and in 3D after ¢ = 2.14. The vertical dashed line in Fig- 
ure 3 indicates the moment of initiation of the 3D simulations. 
Cycle accretion in the 2D simulations begins at ¢ ~ 1.4 and is 
clearly seen as a sequence of spikes in the time dependence of 
Min (see Fig. 3a). Spikes, which are related to the same cycle ac- 
cretion, are also observed in the variation of magnetic flux inside 
R = 2R, (see Fig. 3b). The magnetic fluxes inside the other se- 
lected radii are gradually increased with time because of the in- 
ward advection of the vertical field. The time dependence of 
these fluxes is not significantly influenced by cycle accretion. 

The structure and dynamics of the inner region in model B are 
drastically changed in the 3D simulations. Shortly after the initia- 
tion of the 3D simulations at t = 2.14, the axisymmetric distribu- 
tion of mass near Rm undergoes the Rayleigh-Taylor and, possibly, 
Kelvin-Helmholtz instability (see Kaisig et al. 1992; Spruit et al. 
1995; Chandran 2001; Li & Narayan 2004) with the fastest growing 
azimuthal mode number m ~ 50 (the latter is probably determined 
by our grid resolution). As a result, the empty region inside R, is 
quickly filled, on the free-fall timescale estimated at R,,, with a large 
number of density spikes that move almost radially toward the 
center. These spikes quickly disappear into the black hole, and at a 
later time the inner disk structure is modified toward establishing a 
dynamic quasi-steady state. This state is characterized by a low 
m-mode (m ~ 1—5) spiral flow structure that results from the in- 
teraction of the nearly circular rotating mass of the outer disk with 
the strong vertical magnetic field at Rm. 

Note that a similar low m-mode flow structure was found in the 
simulations of accretion flows onto a magnetic dipole (Romanova 
& Lovelace 2006). The dominance of the low m-modes in the 
quasi-steady state can be explained by the formation of Rossby 
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Fic. 1.— Distributions of density in the meridional plane in model B at two successive moments, (a) t = 5.1153 and (b) t = 5.1458 from axisymmetric 2D simulations. 
The black hole is located on the left, and the small open circle there corresponds to the inner boundary around the black hole at Rin = 2R,. The axis of rotation is in the vertical 
direction. The domain shown has a radial extent 100R, along the equatorial plane and represents a fraction of the full computational domain with Rout = 220R,. The time is in 
units of the orbital period at Rin; = 210R,. The color bars on the bottom indicate the scales for log p (in arbitrary units). The model develops a thick turbulent accretion disk, 
which is seen as the nonuniform mass concentration near the equator. The low-density polar regions are filled with a strong vertical magnetic field. This field causes cycle 
accretion into the black hole. (a) Period of halted accretion, when the mass is accumulated at the magnetospheric radius R,, ~ 15R,. (b) Accretion period. 


vortices near the inner edge (~R,,) of the outer accretion disk 
(see Tagger & Pellat 1999). 

The nonaxisymmetric inner flow is highly time variable and 
experiences quasi-periodic behavior. Figure 4 shows an example 
of the flow structure inside 50R, in model B at two successive 
moments, t = 2.2767 and 2.2867. The flow is essentially 3D 
inside the magnetically dominated region limited by the radius 
~35R, and remains almost axisymmetric on the outside of this 
radius. In the magnetically dominated region, the flow forms mod- 
erately tightened spirals of dense matter, which are clearly seen in 
Figures 4a and 4b. This matter is quickly accreted into the black 


hole with the radial velocity, which is approximately half of the 
free-fall velocity. Such a relatively fast infall is explained by the 
efficient loss of angular momentum by the mass during its inter- 
action with the vertical field. As a result, the accreting mass has a 
sufficiently sub-Keplerian angular velocity, approximately half of 
the Keplerian one. The field is distributed nonuniformly in the disk 
plane, concentrated in bundles that penetrate through the plane in 
very low density, magnetically dominated (with 6 ~ 0.01) regions, 
or magnetic “islands.” The rotating mass interacts with the mag- 
netic bundles and forces them to twist around the disk’s rotational 
axis. In the simulations, this twist is observed as the rotation of 


(b) 


Fic. 2.— Snapshot of magnetic lines in model B at the same moments and on the same spatial scales as in Fig. 1. The component of the field lines parallel to the 
meridional plane is shown. The accretion disk transports the vertical magnetic flux inward, which accumulates in the vicinity of the black hole. Small-scale magnetic loops 
are the result of turbulent motions in the disk and disk corona. (a) Period of halted accretion (see Fig. 1a), in which most of the accumulated magnetic flux is outside the 
black hole horizon. (5) Accretion period (see Fig. 15), in which all the accumulated flux goes through the horizon. The lines have been plotted using the method of Cabral & 


Leedom (1993). 
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Fic. 3.— Evolution of mass accretion rate and magnetic fluxes in model B. 
The time is in units of the orbital period at Rinj. The simulations are performed in 2D, 
assuming axisymmetry, from f = 0 to 2.14 (the latter moment is indicated by the 
vertical dashed line), and in 3D after that. (a) Accretion into the black hole begins 
att ~ 1.3. Starting from t ~ 1.4, the cycle accretion develops (seen as a sequence 
of spikes). In 3D, the cycle accretion disappears. (b) Magnetic fluxes (in arbitrary 
units) through the equatorial plane inside five fixed radii: 2R, (=Rin), 25R,, 50R,, 
100R,, and 210R, (=Rinj). The variability of the flux at Rin is related to the 
variability of the accretion rate in (a). 


magnetic islands around the disk axis. The rotational velocity of the 
islands is typically reduced by a factor of ~0.5—0.9 in comparison 
with the velocity of the surrounding accretion matter. This can be 
explained by the resistance of the large-scale vertical field to such a 
twist. The faster rotation of accretion matter and the slower rotation 
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of the magnetic islands produces a shear flow. The shear flow plays 
two roles in our simulations. First, it provides an exchange of 
momentum and energy between the accreting mass and the vertical 
field. Second, the shear flow results in an ablation of the islands 
caused by magnetic diffusivity, making each island a temporal 
structure. An example of the evolution of magnetic islands can 
be seen in Figure 4: the magnetic islands observed as low-density 
spiral arms above the center in Figure 4a are observed below the 
center in Figure 4 after about half a revolution in the clockwise 
direction. In the latter figure, the islands are apparently reduced 
in size due to the ablation. 

The vertical field ablated from the magnetic islands is carried 
inward by the accretion flow and accumulates on the black hole 
horizon. This accumulation results in quasi-periodic eruptions 
of the field outward from the horizon as soon as the field pressure 
overcomes the dynamic pressure of the accreting mass. The erup- 
tions typically take the form of high-velocity narrow streams (in 
the equatorial cross section) of a low-density, magnetically dom- 
inated medium fountained outward from the black hole. In Figure 45, 
four magnetic islands observed as low-density regions inside 
R z 15R, result from such eruptions, and the eruption of one of 
these islands (to the right of center; see also the steam that pro- 
duced 1t) 1s still continuing at the moment shown. In the conse- 
quent evolution, these islands are pushed outward and stretched in 
the azimuthal direction by the accretion flow, and take a spiral 
shape similar to that shown in Figure 4a. 

Model A evolves faster and accumulates a larger magnetic 
flux at the center than does model B. Figure 5 shows the evolu- 
tion of the accretion rate M;, and magnetic fluxes inside different 
radii in model A. Qualitatively, the evolution of these quantities 
is similar to the evolution of those in model B (see Fig. 3). Quan- 
titatively, however, model A demonstrates significantly larger 
time-averaged accretion rates (by about 2 orders of magnitude) 
and longer quasiperiods of cycle accretion (represented by the in- 
tervals between spikes in the time dependence of Min in Fig. 5a) in 
the 2D simulations. By the end of the 2D simulations at f ~ 1.7, 
this model has the maximum R,, ~ (30—40)R,. 


Fic. 4.— Distributions of density in the equatorial plane in model B at two successive moments, (a) t = 2.2767 and (b) t = 2.2867 from 3D simulations. The black hole 
is located in the center, and the central black circle corresponds to the inner boundary at Rin = 2R,. The shown domain has an extent of 100R, x 100R,. The time is in units 
of the orbital period at Rinj = 210R,. The color bars on the bottom indicate the scales for log p (in arbitrary units). The accretion flow rotates in the clockwise direction and 
forms a spiral structure inside R ~ 35R,. The low-density regions correspond to the magnetically dominated islands, through which the vertical field penetrates the disk. 
(a) Developed structure of spiral density inflows and magnetic islands stretched in the azimuthal direction. (b) The moment after about half a revolution of the spiral 
magnetic pattern seen in (a); new magnetic islands inside R ~ 15A, are the result of an eruption of the vertical magnetic flux from the inner boundary. 
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Fic. 5.— Evolution of mass accretion rate and magnetic fluxes in model A. 
The notations used are the same as in Fig. 3. Note the earlier beginning of accre- 
tion, significantly larger accretion rate and magnetic fluxes, and longer quasiperiods 
between the accretion events in the 2D simulations compared to those in model B 
(see Fig. 3). These are the result of a stronger field injection in model A. 


In the 3D simulations, model A experiences an initial transient 
period, similar to the period of the development of the Rayleigh- 
Taylor instability in model B (see above), in which the nonaxi- 
symmetric, small-scale structures quickly appear and disappear. 
Figure 6 shows an example of the developed low m-mode spiral 
structure in model A obtained after the transient period. This struc- 
ture is clearly dominated by the m = | mode. The magnetically 
dominated region is extended up to R ~ 70R,. Note that the 
spiral-density arms seen in Figure 6 are more open than the arms 
in model B (see Fig. 4a). This could be due to the stronger central 
field in model A. 

Model C is our slowest evolving model and, accordingly, shows 
the slowest rate of accumulation of the central vertical field. This 
model has been calculated only in 2D and demonstrates a qual- 
itative similarity to the axisymmetric evolution of models A and 
B. Cycle accretion, which is caused by the accumulated field, 
begins at t ~ 4.8 in model C. The model demonstrates more ef- 
ficient turbulent motions in the accretion flow. This can be attri- 
buted to weaker magnetic fields, which suppress the MRI and 
convection motions to a lesser extent. At the end of the simu- 
lation at t ~ 6, the model has the maximum Rm ~ 6R,. 


3.2. Poynting Jets 


The 3D simulations of models A and B show that the vertical 
field penetrating the central magnetically dominated regions in 
MADS is twisted around the axis of rotation by the rotating ac- 
cretion flows. The field twist generates electromagnetic pertur- 
bations that propagate outward and transport released energy in 
the form of a Poynting flux (e.g., Landau & Lifshitz 1987) 


S = = (ex B) x Bl (7) 


The Poynting flux is distributed nonuniformly in the polar angles, 
basically showing two components: a jetlike concentration of 
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Fic. 6.— Distribution of density in the equatorial plane in model A att = 1.87 
from 3D simulations. The notations used are the same as in Fig. 4. The domain 
shown has an extent of 300R, x 300R,. The model is characterized by a relatively 
large magnetic flux accumulated at the center. The magnetically dominated region 
has an outer radius R  70R,. The accretion flow demonstrates the domination of 
one-arm (m = 1-mode) spiral structure in the innermost region. This structure 
rotates in the clockwise direction and is split onto several arms at larger radial 
distances. There are another two high-density streams in the process of devel- 
opment (to the left of the center). 


the flux near the poles and a widespread distribution of the flux 
in the equatorial and midpolar angle directions (from 0 ~ 7/4 to 
~37/4). 
Figures 7 and 8 show example @ distributions of the radial 
Poynting flux (solid lines), 
B? Br 
Sr = R4 7 7 (vB), (8) 
at six different radii, SR}, 10R,, 25R,, 50R,, 100R,, and 220R,, 
which cover most of the radial domain in models B and A, re- 
spectively. The distributions shown are averaged in the azimuthal 
direction and in time over the interval Ar ~ 0.05, using a set of 
data files stored during the simulations. Outside the inner, mag- 
netically dominated region (at R Z 35R, in model B and R Z 70R, 
in model A), the outward (positive) Poynting flux in the equatorial 
and midpolar angle directions is supported mostly by the rotation 
of the outer, almost axisymmetric disks, in which the poloidal 
field component is frozen. Such a flux is present in both the axi- 
symmetric 2D and 3D simulations. The 3D simulations introduce 
new important features in the Poynting flux distribution: an increase 
of the flux in the equatorial and midpolar angle directions inside the 
magnetically dominated region, and in the polar directions out- 
side this region (see Figs. 7 and 8). The equatorial flux inside the 
magnetically dominated region is generated by the twist of the 
vertical field in spiraling nonaxisymmetric accretion flows. In 
model B, this flux, represented by bumps in the 0 distributions, 
gradually deviates from the equatorial plane toward the poles as 
the radial distance increases (see Figs. 7a—7d). At R Z 100R,, the 
flux is collimated into bipolar Poynting jets (see Figs 7e and 7 f ). 
In the case of model A, the process of jet collimation is less evident 
and somewhat different; but still, one can observe the formation 
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Fic. 7.— Distribution of Poynting Sp (solid lines) and kinetic Fg (dashed lines) 
radial fluxes over the 0-direction in model B. The fluxes are averaged in the azi- 
muthal direction and in time over the interval At = 0.05, beginning from t = 2.27, 
and normalized to the accretion power M,,c2. The distributions are shown at the 
radial distances 5R,, 10R,, 25R,, 50R,, 100R,, and 220R, (=Rou) in the panels 
from (a) to ( f ), respectively, and illustrate the collimation of Poynting flux into 
bipolar jets. 


of narrow bipolar Poynting jets starting from R ~ 10R, and fur- 
ther development of these jets at large radii (see Figs. 85—8 f ). 

Figures 7 and 8 show 0 distributions of the kinetic flux (dashed 
lines), 


v? 
Fg = RP, (9) 
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Fic. 8.— Same as Fig. 7, but for model A. The fluxes are averaged in time over 
the interval At = 0.04, beginning from t = 1.78. 


jets (except in the outermost region in model A; see Fig. 8/). 
Accordingly, the jet velocity is mostly sub-Alfvénic. However, 
the value of the kinetic flux is relatively large, and this is some- 
what in contradiction to our expectations that MADs can develop 
Poynting flux-dominated jets. The problem of the excessive 
kinetic flux in our simulations can probably be explained by the 
action of the numerical magnetic diffusivity (see § 4), which re- 
sults in an unphysically large mass load for the Poynting jets and 
the consequent excessive kinetic flux in them. 

The Poynting jets are powered by the released binding energy 
of the accretion mass. To estimate quantitatively the amount of 
energy going into the jets, we calculate the Poynting jet “lumi- 
nosity” Éja as a function of the radius R: 


BalR) = J IOS, dO, (10) 


where the integration is taken over the solid angles Q occupy- 
ing the polar regions with 0 < 7/4 and 0 > 37/4 (excluding the 
boundary polar cones; see 8 2). For comparison, we also calculate 
the total Poynting luminosity Etot, which is defined analogously to 
Eja, but with the integration in equation (10) taken over the whole 
sphere. Figures 9 and 10 show the radial profiles of the normalized 
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Fic. 10.— Same as Fig. 9, but for model A. The luminosities are averaged in 
time, as explained in the caption to Fig. 8. 
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Ej (solid lines) and Ex (dashed lines) in models B and A, re- 
spectively. The jet luminosity Eje weakly depends on the radius at 
RZ SOR, in both models, and equals 1.5% in model B and 
0.5% in model A. Here, we quantify the Poynting luminosity in 
units of accretion power, Minc?. The total luminosity Etot in- 
cludes the flux from the bipolar Poynting jets and wide equatorial 
Poynting outflow. The latter component of Erot exceeds Ej. by a 
factor of ~3 at large radii (see Figs. 9 and 10). This could be a 
consequence of the simulation design employed (see 8 2), in 
which the disk accretion at outer radii is mostly provided by the 
global Maxwell stresses. . 

The smaller value of the final (at large radii) relative Ej; in 
model A, in comparison with that in model B, can be explained 
by the different structures of the inner, magnetically dominated re- 
gions in these models. Model A has less tightened spiral density 
arms (see Fig. 6), in which the mass accretes with larger radial ve- 
locity, and therefore delivers less energy to the field. 

Note also that the value of Éja in model A takes a relatively 
large finite value right at the inner boundary Ri, (see Fig. 10), 
whereas Eis in model B begins from a small value at Rin and 
gradually increases outward (see Fig. 9). This difference in the 
behavior of Eje can be attributed to the discussed difference of 
the innermost structure in the considered models. 


4. DISCUSSION AND CONCLUSIONS 


We have performed a numerical study of the formation and 
evolution of quasistationary MADs, which are characterized by 
a strong vertical magnetic field accumulated at the disk centers. 
We employ a simulation design in which the poloidal magnetic 
field of one sign is permanently injected into the computational 
domain at the outer numerical boundary, and the unipolar vertical 
field is transported inward by the accretion flow. The accumulated 
field has a significant impact on the inner flow structure and dy- 
namics in both 2D and 3D simulations. In the axisymmetric 2D 
simulations, the field pressure can temporarily halt the mass fall- 
ing into the black hole, resulting in cycle accretion, in which the 
longer periods of accumulation of mass at the magnetospheric 
radius R„ are followed by short periods of accretion. The 3D sim- 
ulations, however, show no axisymmetric cycle accretion. Instead, 
the accumulated field causes mass to accrete quasi-regularly in the 
form of nonaxisymmetric spiral streams and blobs. We have dem- 
onstrated that 3D MADs can be efficient sources of collimated, 
bipolar Poynting jets that originate in the vicinity of the central 
black hole. These jets are both produced and powered by the in- 
teraction of the spiral mass inflows, with the central field split into 
separate magnetic bundles. The efficiency of conversion of the ac- 
cretion energy Minc? into the Poynting jet energy is up to 1.5% in 
our simulations. This estimate may not be accurate (we believe 
it is underestimated) because of the use of a pseudo-Newtonian 
approximation (see § 2). A better estimate of the efficiency can be 
obtained using general relativistic MHD simulations. 

We have presented the simulation results from three models of 
radiatively inefficient accretion disks, which differ in the strength 
of the injected field. In accordance with previous studies (e.g., 
Stone & Pringle 2001), the structure of the outer disks in these 
models is determined by the field strength. In model A, which 
has the largest injected field, the MRI is suppressed, and the ac- 
cretion flow is driven by global Maxwell stresses, which trans- 
port the excessive angular momentum outward from the disk 
through the disk corona. In models B and C, which have smaller 
injected fields, MRI and turbulent motions are developed. The 
turbulence in these models is axisymmetric and partially supported 
by the efficient convection motions resulting from the dissipation 
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of toroidal and small-scale poloidal magnetic fields. These mo- 
tions cause an increase of the disk thicknesses in comparison 
with nonturbulent model A. The 3D simulations of models A 
and B have demonstrated that nonaxisymmetric motions are not 
important in the outer parts of the disks, outside the inner, mag- 
netically dominated region (these parts remain almost axisym- 
metric), but are very important inside this region, resulting in the 
development of spiral accretion flows and bipolar Poynting jets. 

Numerical magnetic dissipations and reconnections result in 
magnetic diffusivity, which influences the structure and dynam- 
ics of our models. The spatial scale on which the diffusivity oc- 
curs in our simulations is determined by the grid size, which 
greatly exceeds the scales of various resistive mechanisms (e.g., 
Coulomb collisions, dissipative plasma instabilities) in the rele- 
vant astrophysical conditions. Therefore, the magnetic diffusivity 
is significantly overestimated in our models and results in an ex- 
cessive slippage of an accretion flow through the magnetic field. 
This slippage reduces the ability of the flow to drag the vertical 
field inward; however, the numerical diffusivity is not efficient 
enough to totally prevent field accumulation at the disk center. The 
numerical diffusivity suppresses MRI on the scales of the grid 
size, and therefore prevents the development of turbulence in the 
outer regions of our models, where the grid size is increased. An- 
other effect of the numerical magnetic diffusivity is the enhance- 
ment of the ablation of magnetic islands, which are found in the 
3D simulations (see § 3.1). To test the sensitivity of our models 
to magnetic dissipations, we have performed a 2D simulation of 
a model that is similar to model A, but has twice as many grid 
points in the R- and 6-directions. This simulation has demon- 
strated the qualitative similarity of the axisymmetric evolution 
of the high-resolution model and model A: both models show the 
formation of accretion disks, accumulation of the vertical field in 
the disk centers, and development of cycle accretion. The high- 
resolution model forms a thinner laminar accretion disk. Unfor- 
tunately, a more detailed quantitative comparison of these models 
meets some difficulties because of the different properties of the 
mass and field injection region (see § 2), which are changed with 
the change of the resolution. 

The main results of our study, the formation of MADs and 
Poynting jets, have been obtained under the assumption of ra- 
diatively inefficient flows, but we believe that these results can 
also be applied to radiatively efficient, dense accretion disks (e.g., 
Kato et al. 1998). The formation of MADs should not be affected 
by radiative cooling as soon as the central field satisfies the 
equipartition condition, 


B? GMp 
87T Rs 


(11) 


where p is the mass density in the innermost region. The radiative 
losses result in a higher p, and therefore a larger B 1s necessary to 
obtain radiative MADs. We expect that the spiral structure of the 
inner, magnetically dominated region, which is qualitatively sim- 
ilar to that found in our simulations, can be developed in the case 
of the radiative disks. Poynting jets should be a necessary attribute 
of the radiative MADs as well. 

The formation of MADs in radiative disks can be applied to 
the explanation of observations of the low/hard state in black 
hole binaries (for a review, see Remillard & McClintock 2006). 
Here, we briefly discuss the basic 1dea of this application and 
leave more quantitative considerations for future works. We as- 
sume, for example, the development of the MAD in the radiation- 
pressure-dominated accretion disk at the subcritical regime (see 
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Shakura & Sunyaev 1973). In such a disk, the radiation (diffusion) 
timescale fraa ~ H*orp/c can significantly exceed the Keplerian 
time tk = 27R??/(GM)!” at small values of the a parameter, 
o € 0.1, in view of the relation 


t 
—E 34a, (12) 
trad 


which follows from the Shakura-Sunyaev solution. Here, we 
denote H to be the disk half-thickness and or = 0.4 cm? g^! to 
be the Thomson scattering cross section. As soon as trad >> fy in 
the outer Shakura-Sunyaev disk, faq can also significantly ex- 
ceed the accretion velocity facer ~ fy in the inner spiral flow of 
the MAD, in view of the relation 


facer 
Bee X R, (13) 
trad 


which is satisfied in accretion flows with a scaling law of the 
accretion velocity vacer x R-!? and H œ R. 

Having trad >> taccr, one concludes that the radiation is trapped 
inside the spiral flow on the accretion timescale, and that there- 
fore this flow is radiatively inefficient. From the point of view of 
an observer detecting the softer part of the spectrum of outgoing 
radiation (below Av ~ 10 keV), the MAD will look like a Shakura- 
Sunyaev disk truncated at the inner radius Ry, which coincides 
with the transition radius between the inner, magnetically dom- 
inated region and outer axisymmetric accretion disk. Typically, 
in observations, Ry is in the interval from a few tens to hundreds 
of R, (e.g., in Cyg X-1; see Done & Zycki 1999); this value is 
consistent with that obtained in our models A and B. Note that 
the estimate of the radiation timescale trag (see eq. [12]) can be 
modified in the presence of convection in the vertical direction 
of the radiation-pressure-dominated disks (see Bisnovatyi-Kogan 
& Blinnikov 1977; Agol et al. 2001). This likely leads to thinner 
and denser disks than those of the Shakura-Sunyaev solution, and 
a further increase of taq. In the presence of a large-scale magnetic 
field, the synchrotron processes and resistive heating are added to 
the viscous heat production and could modify the disk structure. 

The observed spectra of black hole binaries in the low/hard 
state are dominated by the hard X-ray component (e.g., Done et al. 
2007), and this can be explained by the radiation from the hot, 
optically thin magnetized medium that surrounds the accreting 
spiral flows in MADs and in which the binding energy of these 
flows is released. The synchrotron radiation from the magnetized 
medium and Poynting jets (see Goldston et al. 2005) can be used 
to explain the observed radio luminosity in the low/hard state; this 
luminosity is believed to be due to steady jets (e.g., Corbel et al. 
2003; Gallo et al. 2003). 

Our simulations assume accretion disks around black holes 
and can be relevant to objects with relativistic jets containing 
accreting stellar-mass black holes (e.g., microquasars and black 
holes resulting from Type Ib/c supernova explosions and mergers 
of two compact objects) and supermassive black holes (in galactic 
centers). However, we believe that our main results can also be 
relevant to nonrelativistic astrophysical objects in which accretion 
disks and jets are observed. These objects include, for example, 
young stellar objects and accreting stars (e.g., white dwarfs) in 
binary systems. Qualitatively, we expect that the structure and 
dynamics of MADs and Poynting jets are similar in both the rela- 
tivistic and nonrelativistic cases. We expect, however, large quan- 
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titative differences between these two cases because of the different 
energy-density scales involved in the regions of jet formation. 
Black holes, which are capable of launching jets almost from the 
event horizon, can produce ultrarelativistic jets (e.g., McKinney 
2006). Jets from nonrelativistic objects, which have a surface 
radius R, >> Rg, are limited by velocities v ~ (GM ARG 2 c. 
For example, the latter formula gives an upper estimate of the jet 
velocity of ~400 km s^! from solar-type stars. 

The problem of inward transport and amplification of the ver- 
tical field in turbulent accretion disks has been intensively dis- 
cussed in past and recent years (see, e.g., Spruit & Uzdensky 
2005). The solution of this problem can help to discriminate 
models of accretion disks that are consistent with observations 
(e.g., Meier & Nakamura 2006; Schild et al. 2006). Our simu- 
lation results show that the vertical field is transported inward 
and amplified independent of the disk structure, be it laminar or 
turbulent. It is worth noting, however, that magnetic fields in our 
models are imposed and relatively large. The assumed strength 
of these fields exceeds the possible strength of the self-sustained 
magnetic fields that could be developed due to MRI (Sano et al. 
2004). Therefore, these results should be considered with some 
caution, because they do not represent the case of weak vertical 
magnetic fields. 

The strong vertical field in the center of accretion disks can be, 
in principle, a relic field that is inherited from the previous evo- 
lution. This field can appear, for example, in a supernova merger 
scenario (or the merger of two magnetized neutron stars or a 
black hole with a magnetized neutron star; e.g., Berger et al. 2005) 
or in the course of the gravitational collapse of an extended 
“proto” object (e.g., a protostellar cloud or a supernova pro- 
genitor), which produces a significantly more compact object 
(e.g., a protostar or a black hole). In the latter case, the proto- 
object may contain some amount of the poloidal field, which will 
be amplified and accumulated at the center during the collapse. 
After the formation of the compact object, the remaining non- 
collapsed mass can still move inward, forming an accretion disk 
and confining the field in the vicinity of the object. Depending on 
the relative strength of this relic field and the mass accretion rate, 
MADs and Poynting jets can be developed. The considered sce- 
nario can be applied to young stellar objects (e.g., T Tauri stars; 
Donati et al. 2007) and the hyperaccretion model for gamma-ray 
bursts (Woosley 1993; Paczynski 1998). 

The spiral flow pattern in MADs rotates with about the same 
angular velocity at all radial distances; i.e., it rotates almost as a 
rigid body. Such a rotation can result in quasi-periodic oscilla- 
tions (QPOs) in the emitted radiation, if the disk’s axis is inclined 
to the line of view of the observer (e.g., Alpar & Shaham 1985; 
Lamb et al. 1985; Strohmayer et al. 1996; Lamb & Miller 2001; 
Titarchuk 2003; Tagger & Varniére 2006). The frequency of 
these QPOs should be related to the rotation of the spiral pattern, 
whose angular velocity is defined by the radius of the magnet- 
ically dominated region and can be a fraction (~0.5—1) of the 
orbital velocity at this radius. More investigations are required to 
make quantitative predictions about QPOs from MADs. 
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